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How to directly observe Landau levels in driven-dissipative strained honeycomb 
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We study the driven-dissipative steady-state of a coherently-driven Bose field in a honeycomb 
lattice geometry. In the presence of a suitable spatial modulation of the hopping amplitudes, a 
valley-dependent artificial magnetic field appears and the low-energy eigenmodes have the form of 
relativistic Landau levels. We show how the main properties of the Landau levels can be extracted 
by observing the peaks in the absorption spectrum of the system and the corresponding spatial 
intensity distribution. Finally, quantitative predictions for realistic lattices based on photonic or 
microwave technologies are discussed. 


I. INTRODUCTION 

Artificially-prepared honeycomb lattices offer insight 
into parameter ranges that are out of reach for natural 
graphene, suggesting new directions in which to push ex¬ 
perimental investigations in materials science [1]. One 
of the advantages of such artificial graphene is that the 
tunnelling of particles between different lattice sites can 
be controlled in an independent and flexible way. 

A beautiful result from the theory of electron motion 
in honeycomb lattices is that the effect of an elastic defor¬ 
mation of the material, e.g. as induced by a mechanical 
strain, can be described in terms of an artificial valley- 
dependent magnetic field acting on the electrons. This 
effect was first predicted in carbon nanotubes 013 ] and 
soon extended to two-dimensional graphene mi- A suit¬ 
ably designed strain can generate a constant magnetic 
field such that the electronic spectrum exhibits quantized 
Landau levels, as experimentally observed in graphene 
nanobubbles m- As a consequence of the underlying 
Dirac dispersion of the electron band structure, the spac¬ 
ing of the Landau levels follows a square-root law starting 
from zero energy for both positive and negative energies, 
with the sign of the artificial magnetic field opposite for 
electrons in the vicinity of the two inequivalent Dirac 
points. 

In the photonic context, this physics was first experi¬ 
mentally investigated in El using a distorted lattice of 
optical waveguides: the lattice deformation modulated 
the hopping amplitude between neighbouring waveg¬ 
uides, which in turn generated an artificial pseudo- 
magnetic field for photons. Given the propagating nature 
of the optical set-up, evidence of the Landau levels could 
only be obtained in an indirect way from the localized 
edge modes, which reside in the energy-gaps between the 
Landau levels. In this experiment, it was not possible to 
obtain detailed information on the microscopic structure 
of the Landau levels themselves. 

In this paper, we propose an alternative photonic set¬ 
up consisting of an array of cavities with a honeycomb 
geometry, as inspired by the experiments of m and [13- 
ng; the former used a laterally patterned semiconduc¬ 
tor microcavity device for infra-red light, while the lat¬ 


ter were based on an array of coupled dielectric cylindri¬ 
cal resonators sandwiched between metallic plates. In 
contrast to the propagating waveguide set-up of El, 
such cavity arrays are intrinsically driven-dissipative sys¬ 
tems, which allow the use of spectroscopic techniques 
to characterize their eigenstates. In fact, a general re¬ 
sult m is that the different eigenmodes of a driven- 
dissipative system naturally appear as peaks in the trans¬ 
mission/absorption spectra under a coherent incident 
field. When the pump frequency is set on resonance 
with a given mode spectrally isolated from its neighbours, 
the intensity profiles in both real- and momentum-space 
faithfully reproduce the mode wavefunction HE). 

Here we show how this general spectroscopic technique 
can be applied to the case of a honeycomb lattice to in¬ 
vestigate the Landau levels appearing in the presence of 
distortion. Instead of a complex trigonal distortion of 
the lattice, as considered, for example, in Ennum, 
we consider a simpler geometry with a single hopping 
amplitude that is linearly varied m- Even though this 
configuration is hard to realise by mechanically straining 
a real graphene sample, for the sake of simplicity we shall 
refer to it as the strained honeycomb lattice. 

The structure of the article is the following: in Sec. |TI| 
we review the tight-binding model of the honeycomb lat¬ 
tice, we introduce a configuration that can generate an 
artificial pseudo-magnetic field and we analytically study 
the resulting Landau levels. In Sec. m the analytically 
calculated Landau levels are compared with numerical 
diagonalization of the tight-binding model. The steady- 
state under a coherent driving is studied in Sec. El clear 
signatures of the Landau levels are highlighted and re¬ 
lated to the analytically calculated mode energies and 
mode wavefunctions. In Sec. El we extend our model to 
include next-nearest-neighbour (NNN) hoppings showing 
that the key features of pseudo-Landau levels remain un¬ 
changed. The experimental feasibility of our proposal 
is discussed in Sec. [Vl] using realistic parameters from 
available experiments. Finally, conclusions are drawn in 
Sec. HD 
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II. THE MODEL 


A. The honeycomb lattice 


In the honeycomb lattice, there are two inequivalent 
sites in each real-space unit cell. We label these sites 
as A and 5, and specify the nearest neighbour vec¬ 
tors as: = (—a,0), R 2 = (a/2,v^a/2) and R 3 = 

(a/2, —\/3a/2), where a is the lattice spacing. Within 
a tight-binding approach, we first restrict our attention 
to nearest-neighbour hoppings of strength ti along each 
lattice vector Ri. Later on in Sec. |V| this formalism will 
be extended to also include next-nearest-neighbour hop¬ 
pings. 

In a perfect honeycomb lattice, the nearest neighbour 
hoppings are spatially independent and equal ti = t 2 = 
ts. In the A and B sublattice basis, the wavefunction 
takes the form of a spinor, that can be interpreted as 
a pseudo-spin. Setting the bare site energy as the zero 
energy, the bulk Hamiltonian in momentum space has 
the form: 


n = 





( 1 ) 


with an off-diagonal matrix term: 

Vik) = h +t 2 +t 3 . (2) 

By diagonalizing the Hamiltonian in Eq. Q, we obtain 
the well-known band structure £{k) of the unstrained 
honeycomb lattice, which consists of two bands, sym¬ 
metrically located with respect to the energy zero [8]. 

The first Brillouin zone can be taken in the form 
of a hexagon and is delimited by the Dirac points K 
and K'. The three equivalent K points are located at 



FIG. 1: The strained honeycomb lattice, with Nx unit cells 
along the x-armchair direction and Ny along the y-bearded 
direction (in the figure Nx = Ny = 5). A unit cell contains 
two inequivalent lattice points A and B, separated by a unit 
lattice-spacing a. The unit cell is labelled with two indexes 
(z, j). The nearest neighbour hopping strength along the vec¬ 
tor Ri is denoted by ti . The gradient of hopping t\ along the 
armchair direction is schematically indicated by the thickness 
of the green line connecting A and B sites: the thicker the 
line, the stronger the coupling. 


^ ^ three 

equivalent K' points are located at A' = ^ ~ 3 ^a ) ’ 

K' = Af these six special points, the two 

bands touch at zero energy with a linear dispersion. 


B. Strained honeycomb lattices 

In real graphene, strain is mechanically generated by 
applying suitably designed external forces to the sample, 
so that atoms are physically pushed closer together or 
pulled further apart in a spatially-dependent way. This 
results in a corresponding modulation of both the sample 
geometry and the hopping amplitudes [8]. Typical ge¬ 
ometries used in graphene experiments involve graphene 
nanobubbles uni or trigonal deformations of the honey¬ 
comb lattice [6]. Such trigonal deformations were imple¬ 
mented in the photonic context in the coupled waveguide 
experiment of m- In the present work, we take advan¬ 
tage of the wider design flexibility of artificial graphene to 
show how much simpler geometries with a unidirectional 
hopping gradient can be used to generate an artificial 
magnetic field. 

To understand how a magnetic field naturally appears 
in a strained honeycomb lattice, one can follow the avail¬ 
able literature [8] and make use of an extended tight- 
binding model with a suitably designed spatial depen¬ 
dence of the hopping amplitudes on the site indices 
(bi). 

As a first step, we consider constant, spatially inde¬ 
pendent ti t 2 ts. We expand the Hamiltonian 
in Eq. Q around the K^K' Dirac points by setting 
k = (^cc, g'y — ^47r/(3\/3a)), where \q] ^ 1/a and the 
valley index ^ = ±1 labels the two inequivalent K^K' 
Dirac points. At first order in g, we get: 

~ - *3) - j (4ii +t2+ h) 

1 . . a/3 X 

+ -(2ti — t2 — ts) + — ts). 

With a simple manipulation, we can re-write Eq. § 
in the form: 

« -ivf){hq:o + eA^) + v^^{hqy + eAy) (4) 
where the Dirac velocity is generally no longer isotropic, 

(4^1 + ^2 + ^3) + - ^3) /gs. 

= lt(*2 + h) - - is) 

and an artificial magnetic vector potential appears of 
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components such that: 

vf^eA^ = {t2 - ta) .gx 

v^j^eAy = {2ti -t2- is) • 

In the ti = t2 = ts case of an unperturbed honeycomb 
lattice, the artificial vector potential is obviously zero and 
the Dirac velocities in both two x, y directions are 
equal to: 


with a gradient oriented in the x direction. The positive 
dimensionless parameter r quantifies the intensity of the 
spatially inhomogeneous strain and the positions Xi G 
(—La,/2, 1 /a,/ 2 ). We assume that the ^2,3 hoppings in the 
other two directions are spatially uniform and equal ^2 = 
ts = t. 

We choose the hopping in the middle of the ribbon to 
be ti (0) = t. As we require the hoppings to have the same 
sign at all positions along x, the condition ti{—Lxl 2 ) > 0 
at one edge imposes an upper bound on the strain, 


= ^atl2h (7) 

For generic, but still spatially uniform hoppings, ti 7^ 
^2 7 ^ ^3, the band dispersion is distorted. Provided that 
the hopping imbalance is not too strong, the Dirac cones 
are still present and the valley index ^ keeps its meaning, 
but the position of the Dirac cones in ^-space is shifted 
by the vector potential A pni mi- At the same time, 
the different values of the Dirac velocities in the x,y di¬ 
rections vf) 7^ and the appearance of an imaginary 
part in these signal that neither the group velocity nor 
the pseudo-spin are parallel to the wavevector q anymore. 
Moreover, this imaginary part can be understood as off- 
diagonal components in the Dirac velocity recast in a 
tensorial form [ 5 ]. 

In order to have a non-zero artificial magnetic field, one 
must allow for a spatial dependence of the hopping ele¬ 
ments on the site indices as sketched in Figjl} 

the position of the (local) Dirac points in k space then 
varies along the lattice in real space, following the spatial 
dependence of the artificial vector potential A. The arti¬ 
ficial magnetic field B (directed along the ^ direction) is 
then naturally defined as the curl of A and is generally 
associated with a small spatial dependence of the Dirac 
velocity [ 5 ]. 

As this strain-induced artificial gauge field does not 
break time reversal symmetry, the vector potential A 
and the magnetic field B have opposite signs in the two 
^ = ±1 valleys: in this sense, it is not a true magnetic 
field, but rather a pseudo-magnetic field. In the next 
section, we shall review how a spatially homogeneous 
pseudo-magnetic field rearranges the continuous conical 
dispersion around the Dirac points into a series of dis¬ 
crete Landau levels. 


C. Landau levels in a strained honeycomb lattice 

We consider a ribbon of artificial graphene, infinite 
along the ^-direction and of finite size Lx along the x- 
direction, with Nx unit cells oriented as in Fig .[2 Along 
X, the ribbon is terminated on both ends with bearded 
edges. We assume that the hopping is positive and 
has the following spatial dependence: 


rLx 

6a 


< 1 . 


(9) 


This condition automatically implies that ti{Lxl 2 ) < 2 t 
at the other edge, which guarantees (within a local band 
picture) that no local Lifshitz transition to a gapped 
state [a [22] takes place in the considered ribbon. 

The form of the strain in Eq. (§ leads to a vector 
potential A in the Landau gauge with a constant valley- 
dependent pseudo-magnetic field B = dxAy — dyAx of 
strength: 




The corresponding magnetic length is: 


Ib = \/h/\eB\ = Zajy/^. 


( 10 ) 


( 11 ) 


Inserting the values of the electron charge and the lat¬ 
tice spacing of real graphene, the field in Eq. ( 10 ) would 
correspond to a magnetic field B ^ 2 rl 0 ^ T. 

Inserting the specific form of strain Eq. ^ into Eq. ^ , 
the two components of the Dirac velocity become: 


~ + rtx/3/i 


Vd = 


( 12 ) 


If the modulation of the hopping along the x direction 
is small compared to the hoppings in the other two di¬ 
rections, we can begin our discussion by neglecting the 
spatial variation of the Dirac velocity. In more quanti¬ 
tative terms, the variation of vf) on a magnetic length 
Ib remains small compared to vd as long as the strain 
is weak ^/r ^ 1 . Note that in the strained geometry 
considered in this paper the spatial dependence (12) of 


the Dirac velocity has no impact on the vector potential 
amplitude Ax,y given in Eq. 

The procedure to diagonalize the Hamiltonian is 
closely analogous to the usual one for charged massive 
particles in a uniform magnetic field in free space [8] . At 
the lowest order in q and x, we can express the spatial 
dependence of the hopping amplitudes in terms of a po¬ 
sition operator x and the off-diagonal matrix element in 
Eq. © becomes an operator. 


V = —ihqxVD + fkqyVD + tTxjZa 


(13) 


=tl{Xi) =t (1 + ^ 'r) 


acting on the spatial wavefunction. Using the standard 
(^) canonical relation [x, fx] = i, it is straightforward to show 
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that the commutator [V^, V] =t^T. We can therefore in¬ 
troduce a creation operator d) = V j(t^) which satisfies 
the commutation rules [a, a'*'] = 1 of a quantum harmonic 
oscillator and recast Eq. ([T^ into the form of the creation 
operator of a shifted harmonic oscillator: 


^(\/W<* - *»> “ <“> 

where m = t^/r/(2v]j\ the oscillation center is shifted in 
the X direction by the g^-dependent distance 

ZihVDd ^ 2 ^ /i.N 

xo =- ——qy = -C Ib % (15) 

and the oscillator frequency 

^ ^ ( 16 ) 

recovers the cyclotron frequency of a massless Dirac par¬ 
ticle in a magnetic field, 


: vd' j2\eB\/h. 


(17) 


To obtain the eigenfunctions (j)A^B{x^y)^ we now re¬ 
write the Hamiltonian in Eq. 0 around the Dirac points 
as: 


n 


0 A 

(pB 


= hw 


0 a 
at 0 


4 >a 

(l>B 


= £ 


(f>A 

4>B 


Separating the equations, we have: 

{huj)‘^aa^ (pA = £^(pA 

(hjjj)‘^d) apB = £^(pB- 


(18) 


(19) 

( 20 ) 


Erom Eq. (20), we immediately see that pB is an eigenvec¬ 
tor of the number operator: d)a\N) = N\N) with an non¬ 
negative integer eigenvalue > 0. Therefore, is a ID 
harmonic oscillator eigenfunction \N) with frequency cj, 
centred at the g^-dependent position ( [T^ . Most remark¬ 
ably, for each pB with a given N 0, two independent 
eigenstates exist with opposite energies £ = ±Jiw^/N. 

Obtaining the corresponding pA requires a bit more 
care: from Eq. (18), for N = 0 one finds a single eigen¬ 


state with = 0, while 


Vn\N-1) = ±|Ar- 1) (21) 


Eor each eigenstate, the total wave function in the posi¬ 
tion representation is: 

'i/i-ui \(x v) = “ 1)^ 

^±\n\\ ^y) yp^lx^j Y (^ 11 ^ 1 ) / 

(23) 

where |A^) can be explicitly written in the position rep¬ 
resentation as a Hermite polynomial of degree > 0, 


{x\N) (X H 


N 


X — Xo 

Ib 


(24) 


and we have implicitly assumed that | — 1) = 0. This full 
wavefunction pj{x, y) is therefore a spinor with a Landau 
level wave function in each component: for \n\ > 0, the 
relative sign of the two components pA and 0 b is opposite 
for opposite eigenstates at ±n. Eor n = 0 the state is 
completely localized on the B sublattice. 

Notice that there is no dependence on the valley index 
^ in the spinor of Eq. (23); therefore the wave function 


in the strained lattice is the same for the two Dirac val¬ 
leys. This is an important difference to the case of a 
non-strained system in a (real) external magnetic field, 
where the wave function is valley-dependent and the role 
of the A and B sublattices is inverted when passing from 

? = 1 to ^ =-1 m. 

All the discussion so far neglects the dependence of 
the Dirac velocity on the hopping amplitudes in Eq. 0- 
The first correction to the relativistic Landau levels in 
Eq. ( [^ comes from the spatial dependence of the Dirac 
velocity in Eq. (12). Substituting the value 


trxojSh evaluated at the oscillation center xq into (13), 
one obtains an expression of the Landau levels that in¬ 
cludes the first order correction, 


£n = ±t^/r\n\ y/l- ^Qya. (25) 

This correction shifts each level around the K, K' points 
with a square-root dependence on qy\ the fact that the 
resulting levels are no longer flat is a key difference with 
respect to the standard Landau levels in the presence of 
a real magnetic field, where the Dirac velocity remains 
independent of position. Eurther corrections coming from 
second and higher order terms in the expansion 0 of 
in powers of g go beyond the scope of the present 

work. 


III. COMPARISON WITH NUMERICAL 
RESULTS FROM EXACT DIAGONALIZATION 


for A > 0. 

To summarize, both the positive and negative energy 
states can be organized in a single sequence labelled by 
an integer — oo < n < oo and their energies follow the 
square-root law of relativistic Landau levels with a cy¬ 
clotron frequency cj, 

£n = |n|. (22) 


The analytical results reviewed in the previous section 
were based on several approximations, in particular the 
very notion of an artificial vector potential A relied on a 
local band structure for each value of the hopping 
In order to verify the quantitative accuracy of this ap¬ 
proach, we have numerically diagonalized the full tight- 
binding Hamiltonian and compared the outcome with the 
analytical approximation. 
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FIG. 2: The structure of levels around zero energy as a func¬ 
tion of ky, in units of the bare hopping t. This is numerically 
calculated from exact diagonalization of the tight-binding 
Hamiltonian for a ribbon of Nx = 601 unit cells along x with 
a strain of r = 0.005, and periodic boundary conditions along 
the y-direction. Relativistic Landau levels appear around the 
Dirac points K and K' for kya ~ ±1.21 and ±2.42. Red lines 
indicates the analytical prediction for the lowest Landau lev¬ 
els including the first order correction to the Dirac velocity 
according to Eq. (25). 


sites are separated. The left panel shows the states for 
n = 1: according to Eq. ( [^ , the total wave function 
is '0n=i oc (|0), |1)). On A sites, this corresponds to a 
zero-th order Hermite polynomial, so the wave function 
has a Gaussian profile, as visible from the blue curve. 
On B sites, instead, the wave function is proportional to 
a first-order Hermite polynomial, and has one node as 
visible from the red curve. In the right panel of Fig. 
which shows the n = 2 Landau level, there is one node 
in the ^d-site wave function as expected, and two nodes 
in the 5-site state. The slight asymmetry around the 
centre depends on the sign of the hopping gradient. 

As predicted by the analytical model, exactly at the 
iF, K' points the wave function is symmetrically centred 
in the ribbon while, as ky moves away from the Dirac 
point, its center xq is shifted according to Eq. ( p!^ . The 
finite extension in ky of the Landau level structure around 
the Dirac points, roughly indicated in Fig.j^by the exten¬ 
sion of the red lines, is limited by the Landau level touch¬ 
ing the physical boundary of the system at x = ±I/a,/2. 
In the spectra, this is apparent as a sudden increase of 
the level energy. 


An example of such a comparison is shown in Fig. 
for the energy spectra of a large system of Nx = 601 unit 
cells along x with a relatively small strain parameter of 
r = 0.005. As we are considering periodic boundary 
conditions along y, we can plot the spectrum as a func¬ 
tion of the corresponding ky. The energy levels at zero 
energy are doubly degenerate, consisting of the n = 0 
Landau level and of the localized edge state associated 
with the bearded edges [23]. Around the Dirac points, 
we highlight the formation of quantized Landau levels. 
Their energies are in good agreement with the analyti¬ 
cal prediction of Eq. (25), including the corrections due 
to the spatial dependence of the Dirac velocity The 
slight discrepancies that are visible to a careful eye can 
be explained by the approximations in our analytical cal¬ 
culations, e.g. the neglect of higher-order terms in the 
expansion 


As expected, the agreement between the analytical 
model and the full numerics gets better for smaller val¬ 
ues of the strain parameter r: in this regime, the mag¬ 
netic length extends for a larger number of sites (propor¬ 
tional to 1 /a/t) and the /c-space wavefunction gets more 
localized in the vicinity of the Dirac point. This makes 
the continuum approximation underlying the analytical 
model more and more accurate. At the same time, the 
characteristic value of the vector potential A within the 
real-space wavefunction decreases as y/r, thus reducing 
the importance of the corrections to the isotropic conical 
Dirac dispersion. Of course, this accuracy comes at the 
price of a reduced energy spacing of the Landau levels, 
again proportional to y/r. 

The spatial structure of the eigenfunctions is studied 
in Fig. where the square modulus of the numerical 
eigenfunctions is plotted for two Landau levels with n = 1 
and n = 2, and the contributions from the A and B 


IV. STEADY STATE IN A COHERENTLY 
DRIVEN DISSIPATIVE LATTICE 

After reviewing the main properties of Landau levels 
in a strained honeycomb lattice, we can now describe 
our proposal to probe the Landau spectra by studying 
the steady state of a driven-dissipative system, such as 
a photonic lattice made of a coupled cavity array [12] or 
microwave resonators [T3HIS]. In both cases, the nearest 
neighbour hopping is due to the spatial overlap between 
modes localized on adjacent sites, so the required spatial 
dependence of the hopping Eq. (§ can be tailored by a 
careful design of the distance between neighbouring sites. 
Note that obtaining a gradient of ti along the x direction 



EIG. 3: Square modulus of the wavefunction of Landau lev¬ 
els around the K' point as numerically obtained from exact 
diagonalization of the tight-binding Hamiltonian in Eq. 0, 
choosing ky = 47r/(3A/3tt). The wave function has been sep¬ 
arated for the two sublattices |0 a|^ (blue) and (red). 

Left panel: n = 1, right panel: n — 2. Dots correspond to the 
numerical eigenfunction, lines to the analytical wavefunction 
in Eq. p4|. 
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does not involve distorting the lattice, but only varying 
the distance along neighbouring sites along x. 

Damping naturally appears due to the unavoidable 
photon absorption and radiative losses, so a coherent 
driving by an external coherent light source is considered. 
A related approach was used in [24] to study anomalous 
and integer quantum Hall effects in a square lattice in the 
presence of a synthetic magnetic field and the anomalous 
quantum Hall effect in a honeycomb lattice with sub¬ 
lattice asymmetry. 

We consider a large but finite honeycomb lattice, as 
sketched in Fig. We assume that photons are lost uni¬ 
formly from all sites at the rate 7 and the coherent pump 
is spatially localized in the central part of the lattice. 
This is beneficial to focus on the bulk properties of the 
lattice and to suppress spurious effects due to refiections 
from the lattice edges. 

We coherently pump the system with a monochromatic 
field at frequency ujq. The pump has a spatial amplitude 
fij so that, in the steady state, the time dependent fields 
over A OT B sites at position are: 



FIG. 4: Dots show numerically calculated spectra of the total 
field intensity in Eq. ( |29[ ) as a function of the pumping fre¬ 
quency ujQ for two values of the strain parameter, in arbitrary 
units. The loss rate is h'y/t = 0.005, and we take A = 601 
and ax/a = ay/a = 10 for all panels. Panel A is for r = 0 
(no strain), while panel B is for r = 0.005 (strained) in dark 
blue and r = 0 in light blue. The solid lines are a guide to the 
eye and the red dashed lines indicate the analytical energies 
of the first six Landau levels as predicted by Eq. (22). 


In Fig. [^ we show examples of spectra of the total in¬ 
tensity summed over all lattice sites 


a,,,(T) = aij hj{T) = b^j . ( 26 ) 

with the time-independent amplitudes aij and bij satis¬ 
fying a linear system of Heisenberg equations m- 


h (cJq + 17 ) t bi-ij-i + t bi-ij-^i — fij 

h (cjo + iy) bij -h t\ aij + t + t = fij 

(27) 

for the indexing given in Fig.[2 The strong dependence of 
the non-equilibrium steady-state on the pump frequency 
CJO: as visible in (27), is the key feature allowing for 


the spectroscopic study of the eigenmodes. Conversely, 
waveguide experiments, such as the ones in m, which 
are based on a conservative propagation equation, can 
hardly resolve the different eigenmodes. 

As we want to probe Landau levels arising from strain, 
it is important to separate the contribution of the two 
inequivalent points K and K' as around these points the 
pseudo-magnetic field has opposite signs. To this end, 
we adopt the technique proposed in [24| and we assume 
the coherent driving to have a Gaussian spatial profile 
of width in the two x, y directions, and to have an 
in-plane wave vector in the vicinity of, e.g., the K' Dirac 
point: 


= ( 28 ) 

where Rij is the position vector of the appropriate site 
of the hexagonal lattice and the origin is assumed to be 
located in the middle of the central unit cell. Provided 
the spatial extensions ax^y ^ a, the coherent pump selec¬ 
tively addresses a small region in ^-space in the vicinity 
of the desired K' Dirac point and efficiently excludes the 
other Dirac point. 


It = '^{\ai,j\'^ + \bij\^) (29) 

ij 

as a function of the pump frequency ujq. Depending 
on the specific configuration HMsi, such total inten¬ 
sity spectra can be experimentally accessed with either 
an absorption or a transmission measurement. 

Panel |4]4 shows the spectra of the unstrained case 
r = 0. In this case, the eigenstates of the honeycomb 
lattice [ 8 ] form a continuum with a conical Dirac disper¬ 
sion uj the vicinity of the iF, K' points, so the 

spectrum is a featureless continuum. The small oscilla¬ 
tions in the spectra are due to finite size effects stem¬ 
ming from a small but non-zero reflection at the edges 
and these disappear if larger systems are considered. 

When a small strain r = 0.005 is introduced, pro¬ 
nounced peaks emerge in the spectra shown in panel [4^, 
corresponding to the Landau levels. We compare the nu¬ 
merical spectra with the analytical energies of Landau 
levels, as given by Eq. ( [ 2 ^ . We see that the position of 
the peaks in Fig. is in good agreement with the analyt¬ 
ical values. 

As usual, for a pump frequency close to resonance with 
a peak, the field intensity profile follows the wave function 
of the corresponding mode: the two cases of unstrained 
and strained honeycomb lattices will be separately dis¬ 
cussed in the next sections. 


A. Perfect honeycomb lattice 

For the unstrained case, figureshows the field ampli¬ 
tude \(iij\'^ and \bij\^ on the A and B sites in the steady 
state for two different pumping frequencies cjq and for 
two different spatial extents of the pump along y. The 
first row is obtained for ax/a = 10 and Gy/a = 50, at 
cjo = 0 (panel A) and at cjo > 0 (panel B). The second 
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FIG. 5: Fiel d intensity distribution numerically calculated 
from Eq. (27) for different pumping frequencies uq. N = 601, 
hriji — 0.005, T = 0 and cTx/a = 10 for all panels. Panels 
A and B are for cr^/a = 50, while panels C and D are for 
<jyla — 10. The pumping area spot is highlighted by the 
cyan circle. Panels A and C are for huolt = 0, corresponding 
to the Dirac point, while panels B and D are for a positive 
detuning hwolt = 0.1. 


row is obtained for (Tx/cl = = 10, at the same two 

pumping frequencies. 

The intensity patterns shown in the figure display an 
interesting structure that can be qualitatively understood 
as follows. The pump excites waves that expand radially 
in all directions with the Dirac velocity for a distance of 
the order of vdIj, the so-called conical diffraction [25] . 
The angular dependence around the pump spot is deter¬ 
mined by the matching of the phase profile of the pump 
with the relative phase of the A and B sites at different 
wavevectors in the vicinity of the Dirac point. As can be 
seen in panels A and C for cjo = 0, constructive interfer¬ 
ence reinforces the intensity in the positive-^ direction, 
while destructive interference reduces the intensity in the 
negative-^ direction. 

At finite frequency cJo 7^ 0^ another mechanism con¬ 
tributes to the determination of the pattern: as one 
can see in panel B, the intensity is now concentrated 
laterally in the positive and negative x direction and 
there is almost no intensity in the positive-^ direction. 
This fact can be explained in terms of the momentum 
distribution of the incident field, which does not over¬ 
lap with the resonant Dirac wave at a finite wavevector 
ky uq/vd ^ ^/o'y in the positive-^ direction. As ex¬ 
pected, this feature no longer occurs for a smaller (jy for 
which ky uJojvD < and a good overlap is again 

possible also in the positive-^ direction. As a result, the 
angular distribution is in this case again maximal in the 
positive-^ direction, see panel D. 


B. Strained honeycomb lattice 


The spatial structure of Landau levels in a strained 
honeycomb lattice are illustrated in Fig. a relatively 
weak strain of r = 0.005 is considered and the steady- 
states intensity patterns are shown for different pump¬ 
ing frequencies uq. Panels 0-4 are for frequencies of the 
first five Landau levels n = 0 ... 4, according to Eq. ( [2^ . 
Parameters for all panels were chosen to give the best 
representation of the Landau level wave functions with 
a realistic form of the pump: in particular, a relatively 
wide pump spot along y was needed to keep the jitter of 
the guiding center xq smaller than the magnetic length 
Ib so as to avoid blurring of the Landau level intensity 
pattern. 

To better understand the mode structure, the sepa¬ 
rated field intensity distributions on respectively the A 
and the B sites are shown in the panels labelled with 
A or B: as expected, the qualitative structure of the in¬ 
tensity profile follows the shape of the eigenstates shown 
in Fig. In particular, the numbers of nodes and the 
width of the field amplitudes coincide with what is ex¬ 
pected from the analytical wave functions Eq. (23). 


The particular spiral-like shape is due to the interfer¬ 
ence between the resonant Landau level and its neigh¬ 
bours, that are non-resonantly excited. By decreasing the 
loss rate, in the limit of h^/t 0, the resonant mode is 
more and more dominant while the non-resonant contri¬ 
bution is suppressed and the spiral behaviour disappears. 


V. HONEYCOMB LATTICE WITH 
NEXT-NEAREST-NEIGHBOUR HOPPINGS 

In the discussion so far, we have considered a tight- 
binding model of the honeycomb lattice with only 
nearest-neighbour hoppings, but for realistic systems, it 
may be important to go beyond this approximation. In 
real graphene, for example, electrons hop to next-nearest- 
neighbour (NNN) sites with an amplitude t' that has 
been estimated as being on the order of five percent of 
the nearest-neighbour hopping t [8]. Similar estimates 
have also been made for the experimental realizations of 
an artificial honeycomb lattice, such as in [IMS]. We 
also note that the introduction of a spatially-dependent 
nearest-neighbour hoppings in these experiments would 
result in spatially-dependent NNN hoppings. Motivated 
by this, we now generalize the results of the previous 
sections to the case of non-zero NNN hoppings. 

The momentum-space Hamiltonian corresponding to 
Eq. 0, for a bulk system in presence of NNN hopping 
is: 


y = (VXik) V*{k)\ 

\v{k) vBk)J 


Ho + H' 


(30) 


where the anti-diagonal matrix Ho gives the contribution 
only from nearest neighbour hoppings, and the diagonal 




















FIG. 6: Numerically calculated steady-state intensity distribution for different pump frequencies c^o on resonance with the 
different Landau levels at the K' point. We have chosen N = 601, h'y/t = 0.005, r = 0.005, ax/a = 10 and ay/a = 50 for all 
panels. Panels 0-4 show the total intensity for frequencies uq = 0, oj, uj\/3, ijj\/4 corresponding to the first five Landau 

levels. The pumped region is highlighted by the cyan circle. While the left panels show the complete intensity distribution, the 
central 0A-4A and right 0B-4B panels isolate the intensity distribution on A and B sites. 
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FIG. 7: The honeycomb lattice with next-nearest-neighbour 
(NNN) hoppings. We assume that the NNN hoppings in the 
vertical direction, denoted by t', are all equal and independent 
of position. The NNN hoppings along the horizontal direction 
t[ (i) spatially depend on the position index i of the unit cell 
(introduced in Fig. according to (32). 


matrix describes the NNN hoppings, as depicted in 
Fig.0 In the strained honeycomb lattice, the vertical 
NNN hoppings denoted by t' will remain constant, while 
all other NNN hoppings become position dependent, as 
denoted by where i is the position index i of the 

unit cell. Then we have: 

v^{k) = t' +t\(i) 

+ 1 ) 

V^(fc) = t' 

+ t\{i + l) 

(31) 

where Di = {0,\/3a), £>2 = (3a/2, ■\/3a/2), = 

(3a/2,-y3a/2) are the NNN lattice vectors. We now 
also assume that the hopping has the same spatial 
dependence as the nearest-neighbour hoppings in Eq. 
i.e.: 

i;(i) = (l + ^r) (32) 

Expanding to first order in small ^ 1/a around the 
Dirac point, as done previously in Sec. m we get: 


/ 2 \ 

Vsiq) ~ -3t' - t'r (^—a; + 1 j 

where we have assumed that the spatial dependence is 
weak such that t[{i) — t[{i + 1) ^ t'; this condition is 
guaranteed when r ^ 1 as we have assumed through¬ 
out this paper. Taking t' t, we estimate the energy 
correction due to the NNN hoppings using first order per¬ 
turbation theory as 

A£ = {4>a 0 b ) n' = -it' - ^xo ( 34 ) 


~ -3i' - t'r ( —X - 1 



kytt 


FIG. 8: The structure of levels around zero energy as a 
function of ky^ in units of the bare hopping t. This is nu¬ 
merically calculated from exact diagonalization of the tight- 
binding Hamiltonian for a ribbon of Nx = 601 unit cells along 
X with periodic boundary conditions along the y-direction, in¬ 
cluding next-nearest-neighbour hoppings with t' — 0.08t and 
strain r = 0.005. Red lines indicate the analyt ical prediction 
for the lowest Landau levels according to Eq. (35). 


where we have used that x = xo + y that 

<pA and (pB are eigenfunctions of the number operator 
N = dtd. 

Including the NNN hopping correction, the Landau 
level energy spectrum is now f= A£ + £n : 


^NNN ^ ^ t^/r\n\^yl - ^Qy 


(35) 


where we have used Eq. (15) to evaluate the oscillation 


center xq in Eq. (34). We notice that exactly at the 


Dirac point, the effect of the NNN hoppings is only to 
contribute a global energy shift. 

In Eig. we plot the structure of low-energy levels as 
numerically obtained from exact diagonalization of the 
tight-binding Hamiltonian for the same parameters as 
in Eig. supplemented by a non zero NNN hopping 
t' = O.OSt. Numerical data are compared to the ana¬ 
lytical prediction given by Eq. (35): good agreement is 


clearly visible in the regions around the iF, K' points. 
We have also observed that all qualitative features of the 
wavefunction remain unchanged by the addition of the 
NNN hoppings. 


VI. EXPERIMENTAL REMARKS 

Erom the experimental point of view, the requirements 
to clearly observe the patterns discussed in the previous 
section are a relatively small value of the loss rate and 
a relatively large size of the lattice. On the one hand, 
stronger loss rates may hinder a clear identification of 
the Landau level peaks in the spectra of Eig. [^and are 
responsible for strong mixing of neighbouring eigenstates 
in the spatial pattern of Eig. 

On the other hand, in smaller lattices spurious reflec¬ 
tions at the lattice edges may occur as well as significant 
distortions of the mode wavefunctions. While propaga- 
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tion towards the lattice edges can be a problem for prop¬ 
agating Dirac waves in perfect honeycomb lattices, this 
issue is much less severe in the presence of strain as the 
Landau levels are spatially localized. As a result, finite 
size effects are negligible as soon as the size of the lattice 
is larger the magnetic length Ib- 

In particular, we note that two consecutive Landau 
levels can be resolved if the separation between them is 
larger than the linewidth, namely < En — £n-i- This 
means that the n-th gap is resolved if: 

(36) 

In order for the Landau wave function not to be dis¬ 
torted by the edges of the lattices, we need the mag¬ 
netic length to be much smaller than the size of the 
system. The total length of the lattice in Fig. is 
Lx = {3Nx — l)a/2. The condition Ib ^ Lxl‘1^ then, 
implies that: 


r > 


6V2 

- 1 


2 


(37) 


In order for the key features of Landau level wavefunc- 
tion not to be blurred, such as the number of nodes or 
the width of the wave function, one also needs that the 
position xq of the guiding center jitters by less than a 
magnetic length under the uncertainty oi Qy Ly\ re¬ 
markably, this imposes a condition Ib Ly analogous 
to (37). Finally, we need to keep in mind the condition 
already given in Eq. which set an upper bound on 
the strain r in order to have a physical > 0 hopping 
at all points of the lattice: 


r < 


12 


^Nx 


(38) 


For the sake of concreteness, we can discuss these cri¬ 
teria having in mind a realistic experiment using pho¬ 
tonic na or microwave [ISHIS] technology: state-of-the- 
art samples are in both cases restricted to relatively small 
lattices, with a few tens of sites in each direction. From 


Eq. (38), this imposes an upper bound to the strain 


T < 0.08. 

Eor a realistic loss rate ^ 0.05 of current exper¬ 
iments, a value r ~ 0.07 of the strain should however 
allow us to resolve the first two gaps between Landau 
levels. This is illustrated in the upper panel of Eig. 
where we show the total intensity spectra as a function of 
the pump frequency: peaks corresponding to the lowest 
Landau levels are clearly visible with an excellent agree¬ 
ment with the analytical prediction of Eq. ( [^ . We also 
show that the spectra are only slightly affected when a 
position-dependent NNN hopping ( [3^ is included with 
a realistic amplitude t' = 0.08t. 

In the lower panels of Eig. we show the intensity 
distribution for a pump on resonance with the peak cor¬ 
responding to the n = 1 Landau level. Independently of 



FIG. 9: Numerical calculations using realistic parameters of 
state-of-the-art experimental setups, namely W = = 51, 

h'y/t = 0.05, T = 0.07, (Jx/a — 5 and <Jy/a — 10. Upper panel: 
spectra of the total field intensity as a function of the pumping 
frequency ujq. The blue curve is the spectra obtained with 
only nearest-neighbour hopping, the orange curve includes 
contributions also from the position dependent next-nearest- 
neighbour (NNN) hopping, with t' — 0.08t. The NNN spectra 
has been globally shifted by the factor from Eq. (34) 

in order to directly compare it with the spectra obtained with 
only nearest-neighbour hoppings. Solid lines are a guide to 
the eye and red dashed lines are the analytical energies of 
the first three Landau levels as predicted by Eq. (22). Lower 
panels: steady-state intensity distributions on the A (panel 
lA) and B (panel IB) sites for a pump frequency c^o tuned 
on resonance with the n = 1 Landau level at the K' point, 
in the presence of NNN hopping processes with a realistic 
amplitude t' = 0.08t. 


the presence of a NNN hopping, the peculiar nodal pro¬ 
file of the mode is clearly visible as a central black stripe 
in the B sites intensity pattern shown in panel B. The 
horizontal dark fringes that are visible in the upper and 
lower part of the image are, instead, a spurious effect due 
to reflections on the edges of the lattice. 


VII. CONCLUSIONS 


In this work we have proposed an experiment ally- 
viable spectroscopic method to study Landau levels of 
a Bose field in a coherently-driven dissipative honey¬ 
comb lattice in the presence of an artificial pseudo- 
magnetic field. Eocussing on an inhomogeneous unidi¬ 
rectional strain configuration, we have shown how the 
Landau levels are clearly visible as peaks in the absorp¬ 
tion/transmission spectra of the device. When the pump 
frequency is tuned in the vicinity of a peak, the intensity 
distribution in the lattice provides direct information on 
the microscopic structure of the Landau level wavefunc- 
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tion. 

Experiments along these lines appear feasible with 
state-of-the-art technology using photonic cavity array 
devices in either the infrared or microwave domains. In 
addition to illustrating the physics of Landau levels for 
relativistic massless Dirac particles, such experiments 
would open the way to studying rich wave-propagation 
physics in distorted honeycomb lattices, including Klein 
tunneling [26] and, on more speculative grounds, rela¬ 
tivistic m phenomena. 
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